package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.Date;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

import com.numericalmethod.algoquant.model.ralph2009.Brandt2005Regression;

/**
 * Trial implementation of section 4 experiment with single shock 
 */
public class OneAssetSimulationTrial3 {
	
	// helper functions
	private static double average(double[] a) {
		double result = 0.0;
		for (int i=0;i<a.length;i++) 
			result+=a[i];
		result /= ((double) a.length);
		return result;
	}
	
	private static double stdev(double[] a) {
		double avg = average(a);
		double sum = 0.0;
		for (int i=0;i<a.length;i++) 
			sum+=Math.pow(a[i]-avg, 2.0);
		sum /= ((double) (a.length-1));
		return Math.sqrt(sum);
	}
	
	// CRRA Utility
	private static double utilityFunction(double w, double gamma) {
		
		return Math.pow(w, 1.0-gamma) / (1.0-gamma);
	}
	
	// CRRA differentiated once
	private static double utilityFunctionD1(double w, double gamma) {
		
		return Math.pow(w, -gamma);
	}
	
	// CRRA differentiated twice
	private static double utilityFunctionD2(double w, double gamma) {
		return -gamma * Math.pow(w, -gamma-1.0);
	}
	
	// (TESTING DISPLAY)
	private static void printComparison(double[] nextExcessReturn, double estimatedOptimalWeight, 
			double mu, double sigma, double riskFree, double dt, double gamma) {

		System.out.println("==== COMPARISON ====");

		double averageExcessReturn = average(nextExcessReturn);
		double returnVolatility = stdev(nextExcessReturn);
		System.out.println("Next Step Excess Return = "+averageExcessReturn);
		System.out.println("Next Step Vol = "+returnVolatility);
		System.out.println();
		
		// close form for Mean/Variance portfolio optimization on CRRA Utility 
		// ???CORRECT????
		double closeFormWeight = 1.0/gamma * ((mu-riskFree)/(sigma*sigma));
		
		// (TESTING DISPLAY)
		// the portfolio return mean and vol using our estimated weights
		double ptfMean = averageExcessReturn*estimatedOptimalWeight+(1.0-estimatedOptimalWeight)*riskFree*dt;
		double ptfVol = returnVolatility * estimatedOptimalWeight;
		System.out.println("Using Simulation Estimated Weight: "+estimatedOptimalWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println();

		// the portfolio return mean and vol using the close form weights.
		ptfMean = averageExcessReturn * closeFormWeight+(1.0-closeFormWeight)*riskFree*dt;
		ptfVol = returnVolatility * closeFormWeight;
		System.out.println("Using Mean/Variance Analytic Weight: "+closeFormWeight);
		System.out.println("portfolio mean = "+ptfMean);
		System.out.println("portfolio vol = "+ptfVol);
		System.out.println("portfolio mean/variance = "+(ptfMean/(ptfVol*ptfVol)));
		System.out.println("====================");
		System.out.println();
		
	}
	
	private double[] optimalWeight;
	private double analyticWeight;
	
	public OneAssetSimulationTrial3() {
	}
	
	public double[] getOptimalWeight() {
		return optimalWeight;
	}
	
	public double getAnalyticWeight() {
		return analyticWeight;
	}
	
	
	/**
	 * Solve for optiomal weights
	 * 
	 * @param rand random generator to be used
	 * @param mu annualized return rate
	 * @param sigma annualized volatility
	 * @param riskFree risk free rate
	 * @param pathCounts simulation counts
	 * @param stepsCount the algorithm requires at least 2 steps
	 * @param gamma CRRA utility parameter
	 * @param dt time step size
	 */
	public void solve(final double gamma, final PriceDynamics priceDynamics) {			
		
		// close form solution
		analyticWeight = 1.0/gamma * ((priceDynamics.mu0-priceDynamics.riskFree)/(priceDynamics.sigmaS1*priceDynamics.sigmaS1));
	
		double[][] optimalWeightMatrix = new double[priceDynamics.pathCount][priceDynamics.timePeriod-1];
		optimalWeight = new double [priceDynamics.timePeriod-1];
		double[] optimalWeightVec = new double[priceDynamics.pathCount];
		double[] optimalWeight2 = new double [priceDynamics.timePeriod-1];
		
		double[] aTPlus1Vec = new double[priceDynamics.pathCount];
		double[] bTPlus1Vec = new double[priceDynamics.pathCount];
		double[] vtVec = new double[priceDynamics.pathCount];

		double[] regATPlus1Vec = null;
		double[] regBTPlus1Vec = null;
		double[] regVtVec = null;
		
		//initialization for t=T-1
		RealVector factor1Vec = new ArrayRealVector(priceDynamics.pathCount); 
		factor1Vec.set(1);//set to 1 for t=T-1
		double[] factorVec = factor1Vec.toArray();
		
		// for data-checking purpose in output CSV file only
		//  not used in computation
		double[] averageExcessReturn = new double [priceDynamics.timePeriod-1];
		double[] excessReturnVol = new double [priceDynamics.timePeriod-1];
		

		double[][] excessReturnMatrix = priceDynamics.excessReturnMatrix.getData();
		
		for(int t=priceDynamics.timePeriod-2; t>=0; t--){
			double[] excessReturnVector = priceDynamics.excessReturnMatrix.getColumnVector(t).toArray();
			double[] divYieldVector = priceDynamics.divYieldMatrix.getColumnVector(t+1).toArray();
			
			for(int i=0;i<priceDynamics.pathCount;i++){
				if(t < priceDynamics.timePeriod-2){
					double factor = optimalWeightMatrix[i][t+1] * excessReturnMatrix[i][t+1] + (1+priceDynamics.riskFree*priceDynamics.dt);
					factorVec[i] *= factor;
				}
				
				//aTPlus1Vec[i] = utilityFunctionD1((1+priceDynamics.riskFree*priceDynamics.dt), gamma)*factorVec[i]*excessReturnVector[i]; //Eq.16
				//bTPlus1Vec[i] = utilityFunctionD2((1+priceDynamics.riskFree*priceDynamics.dt), gamma)*(factorVec[i]*factorVec[i])*(excessReturnVector[i]*excessReturnVector[i]);//Eq.17
				aTPlus1Vec[i] = factorVec[i]*excessReturnVector[i]; //Eq.12, redine At+1 = E[phi(Zt+1)*Ret+1] and Bt+1 =E[phi(Zt+1)(Ret+1 * Ret+1)]
				bTPlus1Vec[i] = gamma*(factorVec[i]*factorVec[i])*(excessReturnVector[i]*excessReturnVector[i]);//Eq.12

			}

			//estimate optimal weight per path per time, i.e. x(m,t)
			System.out.println("Fitting At+1 at t="+t + "...");
			regATPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, aTPlus1Vec);
			System.out.println("Fitting Bt+1 at t="+t + "...");
			regBTPlus1Vec = Brandt2005Regression.estimateYValues(excessReturnVector, divYieldVector, bTPlus1Vec);
			for(int i=0;i<priceDynamics.pathCount;i++){
				//optimalWeightMatrix[i][t] = Math.min(Math.max(-1 * regATPlus1Vec[i]/ regBTPlus1Vec[i], 0),1); //shortsell constraint
				optimalWeightMatrix[i][t] = regATPlus1Vec[i]/ regBTPlus1Vec[i]; //shortsell constraint
				optimalWeightVec[i] = optimalWeightMatrix[i][t];
			}			


			//optimalWeight[t] = average(optimalWeightVec); //as suggested in p14 Eq. (21)
			optimalWeight[t] =  average(regATPlus1Vec) / average(regBTPlus1Vec);
			optimalWeight2[t] = average(aTPlus1Vec) / average(bTPlus1Vec); // the original way in simulation 1.
			
			// also save down the average return for benchmarking in CSV file
			averageExcessReturn[t]=average(excessReturnVector);
			excessReturnVol[t]=stdev(excessReturnVector);
			
			//optimalWeight[t]  = -1* average(aTPlus1Vec) / average(bTPlus1Vec);
			// (TESTING DISPLAY)
			System.out.println("At time step "+t+": Estimated Optimal Weight = "+optimalWeight[t] + "(with fitting), "+ optimalWeight2[t] +"(without fitting)" );

			printComparison(excessReturnVector, optimalWeight[t], priceDynamics.mu0, priceDynamics.sigmaS1, priceDynamics.riskFree, priceDynamics.dt, gamma);
			
			//estimate optimal utility per path per time, i.e. V(m,t), as suggested in p14 Eq. (22)
			if(t==0){
				for(int i=0;i<priceDynamics.pathCount;i++){
					vtVec[i] = utilityFunction(factorVec[i], gamma); 
				}
				
				//estimate V(m,e)
				//regVtVec = LSRegression.estimateYValues(excessReturnVector, divYieldVector, vtVec); 
			}
		}
		
		RealMatrix result = new Array2DRowRealMatrix(priceDynamics.timePeriod-1, 4);
		result.setColumn(0, optimalWeight);
		result.setColumn(1, optimalWeight2);
		result.setColumn(2, averageExcessReturn);
		result.setColumn(3, excessReturnVol);

		
		SimpleDateFormat sdf=new SimpleDateFormat("yyyyMMddHHmmss");
		String dateStr=sdf.format(new Date());
		String outputFilename = "C:\\dev\\MAFS\\mafs6010D\\Project\\strategic-asset-allocation\\single-shock-weights"+dateStr+".csv";
		save(outputFilename, result);
		System.out.println("Result saved to "+outputFilename);
	}
	
    public static void save(String fileName, RealMatrix matrix) {
	    try {
	
	        BufferedWriter writer = new BufferedWriter(new FileWriter(fileName));
	        for ( int i = 0; i < matrix.getRowDimension(); i++){  
	        	for(int j=0; j<matrix.getColumnDimension();j++){
	        		writer.write(matrix.getEntry(i, j) + ",");
	        	}
	        	writer.write("\n");
	        }
	        writer.close();
	    } catch(IOException ex) {
	        ex.printStackTrace();
	    }
	}
    
	public static void main(String[] args) throws Exception {
		
		// CRRA utility parameter
		double gamma = 5.0;

		OneAssetSimulationTrial3 t=new OneAssetSimulationTrial3();
		
		//EW with Momentum + MR
		
		double phi = 0.39; // degree of momentumness of the stock price process, range from 0 to 1
		double mu0=1.15/100;
		double mu1=0.012;
		double sigmaS1=6.28/100;
		double alpha = 0.0094;
		double sigmaX1=-8.43/100;
		double sigmaX2=1.52/100;
		

		/*
		// NOTE: all set to zero.
		// ACTUALLY NO LONGER REQUIRED TO SET THIS VAR1 COEFFICIENT
		// The original price dynamics with momentum already implies VAR1 
		double[][] varCoefficient = { 
				{0.0, 0.0, 0.0},		// for log St 
				{0.0, 0.0, 0.0},		// for Mt (not yet supported)
				{0.0, 0.0, 0.0}			// for Xt (not yet supported)
		};*/
		
		
		int shockStartTime=0; // array indexing problem? if start and end = 0, it has no effect....
		int shockEndTime=5;
		// true = use Gaussian random. 
		// 		That will give reasonable variance in simulation, this the weight should be smaller. (i.e. less aggressive)
		// false = use 1.0 in the stock price dZ(t).
		//		that will give near zero variance in simulation and the optimal weight would be very large (since no risk)
		boolean shockWithGaussian=true; 
		
		SingleShockPriceDynamics pd = new SingleShockPriceDynamics(
				mu0, mu1, sigmaS1, alpha, sigmaX1, sigmaX2, phi,
				shockStartTime, shockEndTime, shockWithGaussian);
		
		
		//pd.pathCount=3;
		
		
		pd.simulatePath();
		t.solve(gamma, pd);
		

		System.out.println("Simulation Finished");
	}
	

	


	
	
}
